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We study the dynamics of networks with inhibitory and excitatory leaky-integrate-and-fire neurons 
with short-term synaptic plasticity in the presence of depressive and facilitating mechanisms. The 
dynamics is analyzed by a Heterogeneous Mean-Field approximation, that allows to keep track of the 
effects of structural disorder in the network. We describe the complex behavior of different classes of 
excitatory and inhibitory components, that give rise to a rich dynamical phase-diagram as a function 
of the fraction of inhibitory neurons. By the same mean field approach, we study and solve a global 
inverse problem: reconstructing the degree probability distributions of the inhibitory and excitatory 
components and the fraction of inhibitory neurons from the knowledge of the average synaptic 
activity field. This approach unveils new perspectives in the numerical study of neural network 
dynamics and in the possibility of using these models as testbed for the analysis of experimental 
data. 


I. INTRODUCTION 

Many of the brain activities emerge as the combined ef¬ 
fect of excitatory and inhibitory components associated 
to synaptic plasticity CHS]. In mammalians, the frac¬ 
tion of inhibitory neurons is close to 20-30% [7], and 
it seems plausible that this value has been determined 
by evolutionary constraints, aiming at the effectiveness 
of brain functions. Recently, an explanation has been 
proposed referring to the possibility that such a rate be¬ 
tween inhibitory and excitatory neurons could optimize 
the performances of a neural network j8j. Neurons in 
cortical area can exhibit quite complex scale-free struc¬ 
tures, where inhibitory neurons can play the role of hubs 
that control and moderate the action of the excitatory 
ones [9|. All of these considerations indicate that models 
of neural networks aiming at reproducing a great deal of 
brain functions should take into account the presence of 
both excitatory and inhibitory neurons, organized on a 
suitable network mm- 

The large number of units and the typical high den¬ 
sity of connections in many brain areas suggest that 
a mean field approach can be a proper mathematical 
tool for understanding the large scale dynamics of neu¬ 
ral network models 0 EHnj. Recently, we have ap¬ 
plied a Heterogeneous Mean-Field (HMF) strategy to 
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deal with the dynamics of an excitatory neural network. 
This method retains the basic information on the net¬ 
work topology through the probability distribution P(k) 
of the in-degree density k of synaptic connections at¬ 
tributed to each neuron, and it allows to build the dy¬ 
namics of classes of neurons sharing the same in-degree 
density k, by a suitable discretization of the dynamical 
rule [IB]. The HMF approach is very effective in repro¬ 
ducing the main dynamical features of random dense net¬ 
works of leaky-integrate-and-fire (LIF) excitatory neu¬ 
rons with synaptic short-term plasticity. In particular, 
the structure of quasi-synchronous events and the dis¬ 
tinction between families of locked and unlocked neurons, 
with a rich and complex phenomenology in synchroniza¬ 
tion related to the topological features of the network, 
are fully recovered Un¬ 
interestingly, the HMF approach also allows to solve in 
a natural way a global inverse problem. This consists in 
recovering the unknown degree probability distribution 
P(k) from the knowledge of the average synaptic activ¬ 
ity field. The method has been successfully applied to 
Gaussian and broad degree probability distributions of 
excitatory neurons, and it has been shown to be robust 
with respect to the introduction of noise and disorder 

PS- 

In this paper we show that the HMF strategy can be 
generalized to networks of excitatory and inhibitory neu¬ 
rons, organized in a complex network topology, combin¬ 
ing depressive and facilitating mechanisms. The main 
technical difficulty to overcome is that short-term synap¬ 
tic plasticity obeys different dynamical rules for excita¬ 
tory and inhibitory neurons [Sol EO- Moreover, as dis¬ 
cussed in Sec. [TT[ one has to distinguish between the dy- 




2 


namics of postsynaptic excitatory and inhibitory neurons 
and, for both of these subclasses, between signals coming 
from presynaptic excitatory and inhibitory neurons. 

A comparison between the original network dynamics 
and the corresponding HMF dynamics is reported in Sec. 
III| There we show that, even on a random network with 
both excitatory and inhibitory component, the HMF ap¬ 
proach reproduces the main features of the different com¬ 
plex dynamical regimes. 

In Sec. [TV] we discuss the features of the global synap¬ 
tic activity fields emerging from the HMF dynamics for 
different values of the inhibitory fraction. In particular, 
we show that the system can display a quasiperiodic be¬ 
havior characterized by locked and unlocked neurons, or 
an asynchronous regime where all neurons have differ¬ 
ent oscillation frequencies. Moreover, for a specific value 
of the inhibitory fraction, the system features an opti¬ 
mal synchronization regime, where all neurons display 
the same interspike interval. 

In Sec. [V]we derive an analytic relation between exci¬ 
tatory and inhibitory global synaptic activities, and we 
use this result to show, in Sec. [VlJ that the global in¬ 
verse problem can be solved also for neural networks con¬ 
taining both excitatory and inhibitory components, with 
Gaussian and scale-free degree density distributions. In 
particular, we are able to reconstruct, from the average 
synaptic activity, the degree density distributions and the 
inhibitory fractions on a network with 10% inhibitory 
neurons and two gaussian distributions P{k) for both in¬ 
hibitory and excitatory components. The same holds in 
the case of a network with 30% inhibitory neurons gen¬ 
erated by two P(k) scale-free distributions, with a larger 
average value of k for the inhibitory components. 

Conclusions and perspectives of our research are finally 
presented in Sec. El 

II. LIF EXCITATORY AND INHIBITORY 
NEURONS WITH SYNAPTIC PLASTICITY 

We consider a network of N neurons, either excitatory 
or inhibitory. Calling Vi(t) the membrane potential of 
neuron i. its dynamics is ruled by the LIF model, i.e. 

Ui(t) = a - Vi(t) + I- yn {t) , ( 1 ) 

where a is the common leakage current, and I* yn (t) is 
the synaptic current coming from the connections with 
other neurons. All variables can be rescaled to work with 
adimensional units (see usual). For instance, time is 
rescaled to the membrane time constant r m = 30ms, and 
the spiking threshold Vth of v is set to 1, while its reset 
value is v r = 0. Whenever Vi reaches i>th, the neuron i 
emits a spike and is reset to v r . In our simulations we set 
a = 1.3, so that neurons are in a spiking regime, i.e. even 
in absence of synaptic stimuli they fire periodically with 
a period To = ln(a/(a — 1)). For the coupling dynamics 
we use the Tsodyks, Uziel and Markram (TUM) model, 
a description of short term synaptic plasticity that has 


been successfully tested in experimental setups mm- 
According to [22], the dynamics of the synapse between 
postsynaptic (i.e., receiving) neuron i and presynaptic 
(i.e., transmitting) neuron j is described in terms of the 
fraction of its active, Uij(t), available, Xij(t), and inac¬ 
tive, Zij(t), resources. These quantities are assumed to 
evolve according to the following set of coupled differen¬ 
tial equations: 

Vij(t) = + Uij(t)Xij{t)Sj{t) (2) 

Tin 

Xij{t) = - UijtyXijtySjit) (3) 

.r 

Xij(t) + Vij{t) + Zij(t) = 1, (4) 

where the last equation is a conservation rule and Sj ( t ) = 
JO 8(t — tj(n)) is the spike train of the presynaptic neuron 
j emitting its n-th pulse at time tj (n ). Whenever neuron 
j emits a spike, it activates a fraction u,j of the available 
resources Xi 3 . In between two consecutive spikes, the 
fraction of active resources yi 3 decreases in time with a 
time constant r; n , and the fraction of available resources 
recovers in a time t* the fraction of inactive resources 
Zij. If the postsynaptic neuron i is inhibitory, the re¬ 
covery time is much smaller. In particular, the typical 
phenomenological values are T- m = 0.2, while r® = 3.4, if 
i is inhibitory, and r® = 26.6, if i is excitatory. More¬ 
over, if the index i corresponds to an excitatory neuron, 
Uij(t) is assumed to be constant, namely Uij = U = 0.5, 
otherwise 

= - U P^1 + u (i _ Uii (t))S [I) (5) 

Tf 

where Tf = 33.25 is the facilitation time scale and Uf = 
0.08 is a phenomenological parameter (22] [23] . 

The TUM model equations combine depressive and fa¬ 
cilitating mechanism of plasticity [22] . If the postsynap¬ 
tic neuron is excitatory, the mechanism is purely depres¬ 
sive as a high frequency spiking of presynaptic neurons 
delays the available synaptic resources. If the postsynap¬ 
tic neuron is inhibitory, the dynamics of Ui 3 describes 
a facilitating mechanism, reinforcing the synapse when 
presynaptic neuron j has a high electric activity. Equa¬ 
tions participate in the neural network dynamics 

by specifying in Eq.([T]) the form of the synaptic current 
received by neuron i: 

Ii vn {i) = ( 6 ) 

j#* 

where g is the coupling parameter and the index j labels 
the presynaptic neurons of neuron i. The matrix ele¬ 
ments eij can take the values 0, 1 and -1 if, respectively, 
presynaptic neuron j is disconnected, excitatory or in¬ 
hibitory with respect to postsynaptic neuron i. In our 
first approach, we consider random uncorrelated dense 
networks, i.e. networks where any correlation among dif¬ 
ferent degrees is absent, and the degree is proportional to 
N (this explains the normalization factor 1/iV in Eq.([6]) ). 
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By defining k = k/N the rescaled in-degree, where 
k £ [0, N — 1] is the number of in-connections of a given 
neuron, we can associate to the uncorrelated network its 
in-degree distribution P{k). In general, inhibitory and 
excitatory neurons may have different in-degree distribu¬ 
tions, Pi(k) and P E (k). In particular, Pi(k) and P E (k) 
are the probabilities that an inhibitory and an excitatory 
neuron receives kN inputs from other neurons. In this 
setup, we fix the in-degree distributions of inhibitory and 
excitatory neurons, assuming that typically they have 
randomly distributed outputs. 


III. HETEROGENEOUS MEAN FIELD FOR A 
NETWORK OF EXCITATORY AND 
INHIBITORY NEURONS 

The HMF approach (see [181119) 1 amounts to perform a 
thermodynamic limit, while keeping the neuron in-degree 
density k fixed. Dynamics (JT]) is replaced by an evolution 
rule for classes of neurons labeled by their k 

Vk(t) = a - v~ k (t) + gkY(t) (7) 

where Y ( t ) is the average synaptic activity field. In 
pan! this method was applied to dense uncorrelated 
networks of excitatory LIF neurons, and it revealed a 
very good approximation of dynamics |l]) for any large 
finite network. For instance, the HMF approach is effec¬ 
tive also for sparse uncorrelated networks, provided they 
exhibit a sufficiently large average in-degree m- 

Here we describe how the HMF approach can be ex¬ 
tended to networks made of inhibitory and excitatory 
neurons. In this case, the dynamics of each neuron de¬ 
pends on the number of its inhibitory and excitatory 
presynaptic neurons. This information is stored in the 
adjacency matrix eij, that encodes the network topology. 
Following the HMF strategy, we have to split dynamics 
(|7) into two equations for classes of excitatory ( E) and 
inhibitory (/) postsynaptic neurons with in-degree den¬ 
sity k: 

vf (t) = a - vf(t) + gk(-fiY EI (t) + f E Y EE (t)) (8) 
vI(t) = a-v!(t) + gk(-f I Y II (t) + f E Y IE (t)) . (9) 

The last expressions on the r.h.s. of these equations cor¬ 
respond to the the average synaptic activity fields re¬ 
ceived by postsynaptic E and / neurons with label k 
from their presynaptic E and I neurons. The fractions 
of inhibitory and excitatory neurons are denoted by // 
and f E = l — fi, respectively. 

These equations determine the spike trains .S'- (t) and 
Sj?{t) of inhibitory and excitatory classes of neurons, 

with in-degree density k. On their side, these quanti¬ 
ties enter the set of Eq.s 0, that split into four sets, 
identified by the upperscript (f, *), where both symbols 
can be either I (inhibitory) and E (excitatory), with ” *” 


corresponding to presynaptic neurons with label k. In 
formulae 

Jb *)( t \ 

( 10 ) 

'in 

ij^C t) = - u^\t)x^*\t)Sl(t) (11) 

Tr 

+ + = ( 12 ) 

where u^~’*^=Uifi = E, otherwise 

k 

+ (13) 

Notice that the value of parameter rj depends on the 
type of the postsynaptic neuron. 

These equations can be closed by the consistency rela¬ 
tions defining the average fields that appear in Eq.s ([8]) 
and (|j]): 

p 1 

Y t) * = J P*(k)yf* ] dk . (14) 

For what follows, it is convenient also to define the global 
average fields, Yj and Y E , received by inhibitory and ex¬ 
citatory neurons: 

Yi = —fiYjj(t ) + f E Y IE , (15) 

Y E = —fiY E i(t) + f E Y EE . (16) 


IV. DYNAMICAL EFFECTS OF INHIBITORY 
NEURONS 

In a series of papers □3 DID ED we have analyzed in 
detail the dynamics of random, uncorrelated, dense net¬ 
works of excitatory LIF neurons and successfully com¬ 
pared it with the corresponding HMF dynamics. Anal¬ 
ogously, in this section, we provide a short summary of 
some basic dynamical regimes of the HMF dynamics with 
inhibition presented in Section [TXT] In order to appreciate 
the reliability of this approach, we preliminarily give a 
comparison between the HMF results and direct numeri¬ 
cal simulations, performed on a large finite network made 
of N = 5000 neurons. 

In Fig. [T] we plot the average inter-spike time inter¬ 
val ( ISIf ;), of each neuron as a function of its in-degree 
density k. Data have been obtained for Gaussian proba¬ 
bility distributions Pj(k ) and P E (k ) of k, with (kj) = 0.5, 
(k E ) = 0.7, and standard deviations dj = 0.04 and 
d E = 0.056. Moreover, here and in the following Fig¬ 
ures, we have assumed the phenomenological values of 
the parameters Tf = 33.25 and g = 30 (see [22[ [23)1. 
The matching between the HMF dynamics and direct 
simulations is remarkable. As already observed in fully 
excitatory networks (IBI 119) , excitatory neurons split into 
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FIG. 1: (Color online) Average inter-spike time interval ( ISI ) 
for a network of N = 5000 neurons (violet squares) with 
fi = 0.1 as a function of the in-degree density k. The distri¬ 
butions Pi{k) and Pb(/c) are both Gaussian with (ki) = 0.5, 
(fc B ) = 0.7, 5/ = 0.04 and erg = 0.056, respectively. The 
black (relative to excitatory neurons) and green (relative to 
inhibitory neurons) dots have been obtained by the corre¬ 
sponding HMF dynamics. Notice the plateau region typical 
of the population of excitatory neurons, that is almost absent 
for the population of inhibitory neurons. The inset shows the 
raster plot of the HMF dynamics, where the neuron index s 
is ordered according to the in-degree density k (see text) : 
excitatory neurons correspond to black dots ( 0 < s < 2000) 
while inhibitory neurons to green dots. 


two families, namely periodic (locked) and aperiodic (un¬ 
locked) neurons, respectively observed for k < (fc^) and 
for k > ( kE)■ Inhibitory neurons cover an approximately 
uniform range of higher frequencies. In the inset we re¬ 
port the raster plot (i.e., index of the firing neuron s vs. 
its firing time t) of the HMF dynamics to point out the 
microscopic organization of neurons. Excitatory neurons 
correspond to 0 < s < 2000. In practice, the HMF dy¬ 
namics has been obtained by sampling k with 2000 val¬ 
ues for both groups of neurons, and s has been ordered 
according to decreasing values of k. In fact, the quasi- 
synchronous bursts observed for 0 < s < 1300 are pro¬ 
duced by the locked excitatory neurons in the plateau re¬ 
gion. The unlocked ones (1300 < s < 2000) exhibit quite 
irregular firing behavior, as well as most of the inhibitory 
neurons, that fire more frequently thanks to the facilita¬ 
tion mechanism typical of their synaptic activity. Only 
a small fraction of inhibitory neurons, whose k values 
overlap with those of excitatory neurons in the plateau, 
produce quasi-synchronous bursts, but with a different 
pace and regularity with respect to the excitatory ones. 

We are facing a sort of of higher order locking effect of 
topological origin, induced by excitatory neurons in the 
plateau over equally coupled inhibitory neurons. Despite 
the complexity of the raster plot detailing this regime, 
the average activity fields Ye(I) and Yj(T) exhibit pe¬ 
riodic oscillations, that characterize a dynamical phase 
with a high level of synchrony, mainly driven by locked 
excitatory neurons (e.g., see the inset in the upper panel 
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FIG. 2: (Color online) Maximum (dots) and minimum 

(stars) values of the global field YE{t) as a function of the 
fraction of inhibitory neurons fi obtained from the HMF dy¬ 
namics. The distributions Pis/i(k) and the HMF sampling 
adopted are the same of Fig. [T| 


of Fig. [3]) . 

In what follows, we report how the dynamical regime 
of the HMF dynamics of the model described in Fig. [l] 
changes as a function of the fraction fj of inhibitory neu¬ 
rons. As pointed out in [25], an effective order parameter 
for exploring the HMF phase-diagram is given by the ex¬ 
tremal values of Ypit). In fact, for increasing values of 
//, the amplitude of Yg(t) decreases and the synchrony 
mechanisms inside the network are significantly modi¬ 
fied. In Fig. [ 2 ] we plot the maximum (dots) and min¬ 
imum (stars) values of Yp(t) as a function of //. One 
can distinguish three main regimes. For 0 < // < 0.45 
(regime A), the network dynamics is driven by locked ex¬ 
citatory neurons. The upper panel of Fig. [3] shows the 
average ISI as a function of k for fj = 0.2. We see that 
a large part of excitatory neurons are locked in phase, 
yielding the quasi-synchronous events appearing also in 
the inset of Fig. [T] The inhibitory neurons are mainly 
unlocked and they fire at a higher frequency. In the inset 
we plot the average activity fields received by excitatory 
and inhibitory neurons, Yp{t) and 17(f). Both fields take 
positive values, while 17 (f) > 17 ; (t), as a consequence of 
the facilitation mechanism that increases the synaptic ef¬ 
ficiency during fast firing activity. 

For fj > 0.7 (regime C), Yp(t) becomes negative and 
its amplitude reduces significantly, while both 17 (f) and 
17 (f) do not exhibit any oscillating behavior (not shown). 
This dynamical phase is dominated by inhibitory neurons 
and the natural firing activity of all neurons slows down 
to an irregular behavior, where quasi-synchronous events 
disappear. 

For 0.55 < f 1 < 0.7 (regime B) partial synchroniza¬ 
tion with quasi-synchronous events persist, as in regime 
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FIG. 3: (Color online) Dynamical regimes for different values 
of the fraction of inhibitory neurons // for the same model of 
FigfTl The three panels show the dependence of the average 
ISI j, on k for fj = 0.2 (upper panel), fj = 0.6 (middle panel) 
and fi = 0.5 (lower panel). Green dots stand for inhibitory 
neurons while black dots for excitatory neurons. The insets 
show the average activity fields Ys(t) (black lines) and Yi(t) 
(green dashed). On the left of the lower panel there is an 
additional inset containing the raster plot of firing events. As 
in the inset of Fig. [T] black dots are relative to excitatory 
neurons and green dots to inhibitory neurons. 


A. On the other hand, the microscopic organization of 
firing events is different, as shown in the middle panel 
of FigjSj where fj = 0.6. Looking at the inset, one ob¬ 
serves that inhibitory neurons receive a relatively higher 
amplitude activity field, Y/(f), with respect to excitatory 
ones, Y E (t). The main distinctive feature with respect 
to regime A is that both Y E (t) and Yj(t) take also neg¬ 
ative values, and inhibitory neurons are no more faster 
than the excitatory ones, as one can easily realize looking 
at the average ISI vs. k. In particular, an appreciable 
subset of inhibitory neurons lock at the same frequency 
of locked excitatory ones (see the initial plateau around 
k = 0.4). As to the unlocked inhibitory neurons, they 
fire at lower frequencies with respect to the excitatory 
ones, that, on their side, are locked for large values of k 
(compare the upper and middle panels of Fig§. 

At the edge between regimes A and B, there is a region 
of optimal synchronization (grey band in Figj2]), where 
neither excitatory nor inhibitory neurons prevail. The 
dynamics typical of this region is shown in the lower 
panel of Fig|3j where // = 0.5. The average ISI is in¬ 
dependent of k. This means that all neurons fire with a 
common frequency very close to 1/To, be. the frequency 
of the non-interacting system (<7 = 0 ); only their relative 
phases depend on k . This notwithstanding, the micro¬ 
scopic organization of firing events is still quite complex: 
as shown in the raster plot in the left inset, there is a ma¬ 
jority of inhibitory and excitatory neurons participating 
the same quasi-synchronous events, i.e. they are almost 
in phase. The right inset shows that Y E {t) and Yj(t) ex¬ 
hibit periodic fluctuations of very small amplitude (i.e. 
g ps 0 ), apart narrow activity peaks, that correspond 
to the quasi-synchronous events shown in the left inset. 
This is quite an interesting collective dynamical behav¬ 
ior, emerging in a weakly interacting network, where the 
complex organization of the phases of equally periodic 
excitatory and inhibitory neurons is determined by the 
random structure of the network, i.e. by the in-degree 
density distributions P E (k) and Pi(k). 


V. RELATION BETWEEN AVERAGE 
EXCITATORY AND INHIBITORY ACTIVITY 
FIELDS 


The HMF approach provides also the possibility of 
working out an analytic study of the complex network 
dynamics described in the previous section. For instance, 
direct inspection of Fig. [3] suggests that the average 
synaptic activity fields defined in Eq.s(15l and (16 1 ex¬ 
hibit quite a similar behavior in time, despite single in¬ 
hibitory and excitatory neurons receive different synaptic 
activity fields and emit different spike-trains. In Fig. [4] 
we focus on the time evolution of Y E (t) and Yj(t) for 
the HMF with // = 0.1 (See Figure caption for details). 
Apart from small fluctuations, the two fields coincide by 
a suitable rescaling of their amplitudes (see left inset). 
An analytic estimate of the rescaling factor can be ob- 
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FIG. 4: Time evolution of Yi(t) (green dashed line) and Ye(I) 
(black continuous line) for the HMF dynamics with fi = 0.1. 
In the inset on the left, Tj(t) has been rescaled using the 
common period of the global fields and the factor obtained 
analytically (see text). In the inset on the right we show the 
comparison between the scale factor deduced from simulations 
(black circles) and the one obtained analytically (stars). The 
latter has been calculated by dividing r.h.s. term of Eq. (171 
by r.h.s. term of Eq. (181, where the period T was obtained 
from the time evolution of the global fields. The distributions 
PE/i(k) and the HMF sampling adopted are the same of Fig. 
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Notice that, as the two fields decrease exponentially 
with the same time constant in between two consecutive 
spikes, the field y?’*{t) is equal to i/~’*(t) apart from a 
scaling factor that can be calculated at their maximum 
values, obtained from Eq.s |T7| ) -fll 8 |). In the right in¬ 
set of Fig. [4] we compare the scaling factor computed 
numerically with the analytic prediction obtained from 
Eq.s([T7f ([19]). The agreement is quite good, despite the 
simplifying assumptions introduced in the analytic esti¬ 
mate. Let us point out that, in principle, this should 
hold for a periodic dynamics. Simulations indicate that 
it is effective also when the frequencies of many neurons 
are not too far from the one of the average activity fields 
Fg(t) and Y/(f). 


tained by a heuristic argument. Since Ye and Yj depend 
on the fields and V^’*\ respectively (see Eq.|l4|), 

their difference can be traced bac k to the different dy- 
namical behavior of (see Eq.(13)), that comes into 


play at the firing events. Accordingly, in between two 
spikes, Yj(t) and Ye{£) follow the same dynamics, i.e. an 
exponential decay with the same time constant Ti n . Let 
us consider a neuron that, firing its spike train, gener¬ 
ates postsynaptic fields y'-’* ( t ) and assume that it emits 
spikes at a constant rate, i.e. its synaptic activity field is 
periodic with the same period T of the average activity 
field (actually, locked neurons display such a behavior). 
By imposing the periodicity properties to Eq.s ( |10| ) and 
(131 i.e. yt’*(t) = j/t’*(t + T), we can obtain an explicit 
expression of their time dependence. Both fields exhibit 
the same exponential decay, with time constant Ti n , and 
their amplitudes are found to depend on the different 
boundary conditions at firing events for excitatory and 
inhibitory neurons. In formulae we report their maxi¬ 
mum values, V^'l IAX and yz* MAX ’ achieved immediately 
after the spike emission: 


VI. THE INVERSE PROBLEM WITH 
INHIBITORY NEURONS 

Global synaptic activity fields in extended regions of 
the brain can be more accessible to experimental mea¬ 
surements than single-neuron activities. As shown in 
two previous papers [1510E1 > i n the HMF frame one can 
recover the degree distribution of a fully excitatory LIF 
network from the knowledge of its global synaptic activ¬ 
ity field. In other words, the HMF formulation allows one 
to solve a global inverse problem. Here we discuss how to 
extend such a result to networks made of inhibitory and 
excitatory LIF neurons. The method can be extended to 
networks with different single neuron models. 

Let us assume that we have access to the measure 
of the average activity field received by neurons, i.e. 
Y(t) = f eYe( t) + fiYj(t). We would recover, within a 
reasonable accuracy, P^(fc), Pr(k) and //, i.e. the prob¬ 
ability distributions of the equivalent in-degree density 
k of excitatory and inhibitory neurons, as well as their 
fraction. The procedure goes through the following steps. 

(i) As discussed in SecjVj Yk(£) and Y/(t) can be 
rescaled by a suitable proportionality constant, 7 , whose 
explicit expression depends only on phenomenological pa¬ 
rameters of the model (see Eq.s([T7| (fl9])) . Accordingly, 
we can write 

Y,(t) = 1 Y E (t) = TT rY±L Ii (20) 
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FIG. 5: (Color online) Reconstruction of Ps(fc) and Pi(k) 
for the network of Fig. [T] Continuous curves are the expected 
distributions while red circles (excitatory neurons) and blue 
stars (inhibitory neurons) are the reconstructions obtained 
with the self-consistent inversion equation. 


and we consider Y E (t) and Yj(t) as functions of Y ( t ) and 
of the unknown fraction //; 

(ii) For each value of //, we integrate the dynamics ([8]) 
and (l9|), that produce the spike trains S?(t) that allow 


to integrate Eq.s ( |10| ([13]). 

(iii) We use yto impose the self-consistency 
condition (see (|l4|)) 




( 21 ) 


We can write the equations analogous to (15) and (16) 

Yi = -fiYn{t) + f eYie , (22) 

Ye = ~fiY E i(t) + f eY E e • (23) 

and 

Y(t) = f E Y E (t) + f I Y I (t) 


(iv) Notice that the effective field Y(t) depends on the 
quantities to be recovered, namely P E (k ), Pi(k ) and //, 
and the self-consistency condition © should hold only 
if Y(t) —1- Y{t). In practice, a suitable estimate of the 
unknown quantities can be obtained by minimizing the 
variance 

a 2 = -^— f\Y{t)-Y(t)fdt (24) 
ti ^ t 0 j tQ 

where [to,fi] is the measurement time interval of Y(t). 
The minimization procedure can be achieved by a zero 
temperature Montecarlo algorithm, as for the purely ex¬ 
citatory case (see ,18]). 

The inverse problem procedure (i) (iv) allows to re¬ 
cover quite well the fraction fj = 0.1 of inhibitory neu¬ 
rons. In Fig. [H] we show the reconstruction of P E (k ), 


Pi (k) for the dynamics reported in Fig. [lj This analysis 
confirms that, in the case of Gaussian in-degree density 
distributions for both excitatory and inhibitory neurons, 
the average synaptic activity signal can be efficiently in¬ 
verted. 

Interestingly, this global inverse procedure can be ap¬ 
plied also to the case of broad power law distributions, 
with inhibitory neurons typically displaying higher con¬ 
nectivities jSj- As an example, here we report just the 
case of a network with // = 0.3 and where P E (k ) and 
P/(fc) are power law distributions, scaling as k~ a 0 . In 
order to avoid too small values of k we impose a lower 
cutoff, k ^ and k I ml over both probability distributions, 
that are accordingly normalized. 

In Fig. [6] we report the results of the inverse problem 
procedure. The average synaptic activity held Y (f) has 
been computed from the dynamics of a finite network 
of N = 5000 neurons. In panels A and B we show the 
reconstruction of P E {k) and P/(fc). The minimization 
procedure (step (iv) ) provides quite an accurate recon¬ 
struction of the fraction // ss 0.3 and of the distribution 
P E {k ) over the whole range of definition, while Pi(k) is 
recovered just for k > . This result indicates that a 

more refined algorithm should be employed to improve 
the quality of the inversion. 

We remark that, despite in this case the reconstruc¬ 
tion of Pi(k) is not completely reliable, the presence of a 
fraction of inhibitory neurons in the inversion procedure 
is crucial also for the reconstruction of the excitatory dis¬ 
tribution. Let us consider the global held and implement 
the inverse problem in absence of inhibitory neurons (ac¬ 
tually the procedure reported in [18] ). In the lower panel 
of Fig. © we can observe that, by omitting the presence 
of inhibitory population, the reconstruction of the net¬ 
work structure yields a bad regeneration of P E (k) as well 
(see the inset). Moreover, we compare the reconstruc¬ 
tion performance by plotting the convergence of a in the 
Montecarlo minimization. The curve, in presence of only 
excitatory neurons, converges to higher values of < 7 . More 
precisely, the relative error 6 = a/(Y), which is 48% in 
the case without inhibitory neurons, reduces to 1% when 
inhibitory neurons are taken into account. This implies 
that the matching between Y(t) and Y(t) is neatly im¬ 
proved by considering the inhibitory population. 


VII. CONCLUSIONS AND PERSPECTIVES 

We have studied the dynamics of random uncorrelated 
dense networks of LIF neurons with inhibitory and ex¬ 
citatory components by the Heterogeneous Mean Field 
approximation. This method proves extremely effective 
in reproducing the complex emerging dynamical phases 
of the system and provides significative advantages, both 
in numerical simulations and in the analytic approach to 
the inverse problem, that can be formulated in terms of 
average properties. 









FIG. 6: (Color online) Reconstruction with power law con¬ 
nectivity degree distributions. Panels A and B display the re¬ 
constructions of Ps(fc) and Pi(k ) respectively (black circles). 
The fraction of inhibitory neurons results to be fi = 0.3, as 
wanted. Red continuous line is the power law AT 5 with cutoff 
kin = 0.8 and km = 0.4. The lower panel (where n is the MC 
step) shows the procedure result when the global field is in¬ 
verted by considering only excitatory neurons. In particular, 
we compare the convergence of the parameter a in the case 
of panels A and B where the presence of inhibitory neurons 
is taken into account (black continuous line) with that in the 
case where the inverse problem has been performed by con¬ 
sidering only excitatory population (dashed red line). In the 
inset we show the reconstructed Pg(fc) in this last case and 
we see that it is quite far from the expected one of panel A. 


The model presents a very rich dynamical phase dia¬ 
gram where inhibitory and excitatory components fea¬ 
ture different complex evolutions. Such a complexity 
does not restrain the HMF approach from offering an 
interesting first insight on the global dynamics, that can 
be grasped directly from a simplified version of the HMF 
equations. Precisely, if the dynamics of the two types of 
synapses would be the same, the presence of a certain 
fraction fj of inhibitory neurons could be described by 
introducing an effective in-degree density distribution. 

Indeed, since in this simple case the same activity Helds 
are transmitted to inhibitory or excitatory neurons, we 
just deal with a single set of evolution equations, where 
the dependence on the inhibitory or excitatory nature 
of presynaptic and postsynaptic neurons can be omitted. 
In fact, any neuron with degree k receives a field gkY , 
where 


Y = 


IePe ^) - fiPiik ) yj,(t)dk 


If the term in square brackets has a definite sign, the 
network is equivalent to a completely inhibitory or exci¬ 
tatory (depending on the sign) network with an effective 


probability distribution F(k) 


fEPE(k) - f I P I (k) 


. In 


particular, if Pg = P/ the introduction of a fraction of 
inhibitory neurons fj is equivalent to an effective dilu¬ 
tion in the original network, obtained through 2 fj cuts 
of the links. If the term in square brackets has no defi¬ 
nite sign, P(fc) is not a probability distribution, so that 
the real dynamics does not correspond to an equivalent 
excitatory or inhibitory network. 


In the more complex case considered in this paper, the 
inhibitory dynamics is characterized by a facilitation ef¬ 
fect and the fields received by inhibitory neurons turns 
out to be larger than the excitatory field. This difference 
can be approximately estimated by an analytic argument, 
that allows to implement an approximated inverse prob¬ 
lem for the distributions P*(fc) even in the presence of 
inhibition. 


The global inverse problem proves very effective in re¬ 
producing the in-degree density probability distributions 
of the two populations, together with the fraction of in¬ 
hibitory neurons, for meaningful distributions ranging 
from Gaussians to power laws. On a technical ground, 
the inversion could be improved by adopting more re¬ 
fined minimization procedures, but it is significant that, 
even the simple zero-temperature Monte Carlo method 
adopted here, it is sufficient to provide the wanted out¬ 
come. 


This result paves the way to an analysis of experimen¬ 
tal data of the average synaptic activity fields from broad 
regions in the brain. Moreover, the overall approach can 
be extended to models of correlated dense networks and 
with single neuron dynamics different from LIF. 
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